options.Nbell       = 100;        % Number of Bellman iterations before Newton.
options.Nnewt       = 15;       % Maximum number of Newton steps
options.tolc        = 1e-12;     % Tol on value functions
options.T_irf       = 50;       % Number of periods for IRFs
options.T           = 150;     % Number of periods for MIT shock
options.Terg        = 500;      % Number of periods to burn
options.hpisteps    = 100;

% Stationary
options.L = [];
options.itermaxL    = 5000;     % Max number of iterations to find L
options.tolL        = 1e-12;    % Tol for L
options.tolK        = 1e-5;     % Tol for equilibrium K
options.itermaxK    = 100;      % Max iterations for bisection
options.cresult = [];

options.tolD        = 1e-12;
options.itermaxp    = 25;

options.damp = 0;
options.damp_mit = .5;
options.damp_DC  = .1;

optset('bisect', 'tol', 1e-08)

%% Set globals and parameters
% % try vavra numbers
glob.n          = [50,30];   % Number of nodes for b,a
glob.nf         = [50,30];

glob.spliorder  = [1,1];    % Order of splines (Envelope Condition Method seems only robust when quadratic in k)
glob.Na         = 30;           % Number of nodes for quadrature
glob.curv       = .1;            % Amount of curvature for p grid

glob.Nsignal    = glob.n(2);1e6;

glob.minb = 0;
glob.maxb = 20;


%% Model parameters

% preferences
param.beta      = 0.96;               

%  rho_s    phi_L  sigma epsilon   ce_mu     d_E
% 0.60338 0.076289 2.9199   2.835 -2.3445 0.59226

% idiosyncratic shock
param.rho_s           = 0.79; 
param.sigma_s         = 0.17; 


% labor adjustment cost
param.phi_L           = 0.055;

% Demand parameters

param.omega           = 1;                  
         
param.A = 1;

% entry and exit
param.gamma = .11;     % exogenous exit rate

% fixed cost of production
%param.mu_f     = -3.834;
param.sigma_f  = 1.55887;

% sunk cost of entry 
param.ce       = exp(0); 0.0955;
param.ce_noise = 1e-3;

param.mu_f = -1e6; log(param.ce/2) - param.sigma_f^2/2;

% dist of entrant productivity signals
%param.xi     = .25215;
param.d_E     = .4;

param.nu      = .5; % inverse Frisch elasticity
                     % same as in clementi
                     % can calibrate aggregate shock size to hit moments

% free labor adjustment 
param.delta    = 0.1;

glob.mink       = 0;
glob.maxk       = param.sigma^(param.sigma/param.epsilon);

%% Find initial steady state
options.itermax_DC = 100000;

superelasticities = [.51 .53 .55 .57 .59 .61 .63];
for i = 1:length(superelasticities)

    se = superelasticities(i);
    
    param.sigma           = 20;              % average elasticity
    param.epsilon         = se*param.sigma;     % epsilon/sigma is superelasticity, 0 is CES 

    
    [param,glob]    = setup(param,glob,options);

    eq_SS           = solve_eq_ss(param, glob, options);
    param.M         = eq_SS.M;
    save('eq_SS.mat', 'eq_SS', 'param', 'glob', 'options');
    % param.mu_entry  = eq_SS.mu_entry;
    % param.mu_f      = param.mu_entry - param.mu_f_diff*abs(param.mu_entry);

    glob.agg = kimball_agg(1, eq_SS, param, glob); % Fix this value!
    L        = sum(eq_SS.l.*eq_SS.L);
    param.psi = 1/(L^param.nu); % psi = 1/(L^nu*C)

    % check that C = 1... 
    tosolve = @(C) kimball_agg(C, eq_SS, param, glob) - glob.agg;

    fprintf('C is %.5f \n', bisect(tosolve, .5, 1.5));

    % %% nvm
    % calib = [0.093017 10.753 1.7897 -5.0936  1.1223 2.196 0.62235];
    % mom = compute_moments_for_calibration(calib, param, glob, options);

    %%
    mom = compute_moments(eq_SS, param, glob, options);
    beta_l(i) = mom.beta_l;
end

%%
close all

plot(superelasticities, beta_l,'LineWidth',4)

title('')
xlabel('Super-elasticity')
ylabel('Within Firm Labor-Sales Regression')

set(gcf,'units','points','position',[10,10,1000,600])
set(findall(gcf,'-property','FontSize'),'FontSize',16)

print('-dpng', 'figures/superelasticity_identification.png')
